slope <- map_dfr(rigal_trends,
~.x %>%
select(siteid, linear_slope),
.id = "response"
)
slope_df <- rigal_slp_df
Summary
summary_slope <- slope %>%
group_by(response) %>%
summarise(summ = list(enframe(summary_distribution(linear_slope, na.rm = TRUE)))) %>%
unnest(cols = summ) %>%
pivot_wider(names_from = "name", values_from = "value")
summary_slope %>%
mutate(response = get_var_replacement()[response]) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "bordered", "hover")) %>%
scroll_box(width = "100%", height = "400px")
|
response
|
min
|
1st_quart
|
median
|
2nd_quart
|
max
|
mean
|
sd
|
n
|
n_na
|
frac_na
|
|
Appearance
|
-0.0384615
|
0.0000000
|
0.0053644
|
0.0133333
|
0.0986395
|
0.0084877
|
0.0122101
|
5516
|
0
|
0
|
|
Chao (binary, similarity)
|
-0.1209148
|
-0.0124542
|
0.0000000
|
0.0000000
|
0.0748988
|
-0.0080145
|
0.0165902
|
5516
|
0
|
0
|
|
Chao evenness
|
-1.2609304
|
-0.0247924
|
0.0000000
|
0.0282298
|
2.1442931
|
0.0021943
|
0.0869975
|
5516
|
0
|
0
|
|
Chao species richness
|
-2.8250079
|
-0.0435423
|
0.0027879
|
0.0694623
|
2.2999419
|
0.0160926
|
0.1864729
|
5516
|
0
|
0
|
|
Chao shannon
|
-1.0688101
|
-0.0337486
|
0.0033085
|
0.0500790
|
0.8412016
|
0.0109086
|
0.1104958
|
5516
|
0
|
0
|
|
Chao simpson
|
-2.2048429
|
-0.0294179
|
0.0026013
|
0.0440541
|
0.7745609
|
0.0077291
|
0.0994345
|
5516
|
0
|
0
|
|
Disappearance
|
-0.0300000
|
0.0000000
|
0.0030656
|
0.0097046
|
0.0651261
|
0.0060655
|
0.0093288
|
5516
|
0
|
0
|
|
Evenness
|
-0.1006653
|
-0.0074240
|
0.0000280
|
0.0085548
|
0.1019762
|
0.0010409
|
0.0185918
|
5516
|
0
|
0
|
|
SER_a (rel abundance)
|
-0.1258587
|
-0.0296959
|
-0.0134346
|
-0.0016327
|
0.0634554
|
-0.0184716
|
0.0216062
|
5516
|
0
|
0
|
|
Horn (binary, similarity)
|
-0.1179941
|
-0.0192192
|
-0.0104115
|
-0.0035745
|
0.0421363
|
-0.0131735
|
0.0142384
|
5516
|
0
|
0
|
|
Jaccard (binary, similarity)
|
-0.1209148
|
-0.0268493
|
-0.0152704
|
-0.0059477
|
0.0384615
|
-0.0180388
|
0.0171575
|
5516
|
0
|
0
|
|
Jaccard (binary, dissimilarity)
|
-0.0384615
|
0.0059477
|
0.0152704
|
0.0268493
|
0.1209148
|
0.0180388
|
0.0171575
|
5516
|
0
|
0
|
|
Log species richness
|
-20.6068500
|
-1.0398436
|
0.1305570
|
1.7745696
|
43.8907431
|
0.5690175
|
3.5608475
|
5516
|
0
|
0
|
|
Log total abundance
|
-66.9185877
|
-3.4233489
|
0.5582744
|
4.8750841
|
80.6903559
|
1.2164360
|
8.8216646
|
5516
|
0
|
0
|
|
Nestedness (jaccard)
|
-0.0577731
|
0.0000000
|
0.0053424
|
0.0147321
|
0.1377152
|
0.0090678
|
0.0154355
|
5516
|
0
|
0
|
|
Shannon
|
-0.2007711
|
-0.0106280
|
0.0008848
|
0.0163682
|
0.2212985
|
0.0031355
|
0.0287062
|
5516
|
0
|
0
|
|
Simpson
|
-0.0829149
|
-0.0055055
|
0.0002939
|
0.0073310
|
0.0921231
|
0.0013521
|
0.0142783
|
5516
|
0
|
0
|
|
Species richness
|
-1.6727273
|
-0.0419836
|
0.0038939
|
0.0729567
|
2.0714286
|
0.0233072
|
0.1763604
|
5516
|
0
|
0
|
|
Total turnover (codyn)
|
-0.0384615
|
0.0041152
|
0.0119597
|
0.0220058
|
0.0986395
|
0.0145533
|
0.0147505
|
5516
|
0
|
0
|
|
Total abundance
|
-1726.0511201
|
-1.7400804
|
0.0955316
|
2.3416650
|
1027.3686567
|
2.2395645
|
42.9647716
|
5516
|
0
|
0
|
|
Turnover (jaccard)
|
-0.0769231
|
0.0000000
|
0.0022039
|
0.0156534
|
0.1512605
|
0.0089710
|
0.0161379
|
5516
|
0
|
0
|
bs <- map(var_temporal_trends, ~summary_distribution(slope_df[[.x]]))
names(bs) <- var_temporal_trends
- Jaccard trends: -0.018 (0.017 s.d.)
(-0.01 in @dornelas_assemblage_2014).
- meaning -0.18 % change per decade in average
- Jaccard decreased in 89%
of the sites (79% in @dornelas_assemblage_2014)
- Richness trends:
- +0.6% richness change per year
(3.6 s.d.), i.e.
6% per decade.
Check trends
p <- map2(rigal_trends, names(rigal_trends),
~.x %>%
ggplot(aes(x = linear_slope)) +
geom_histogram() +
labs(title = .y)
)
plot_grid(plotlist = p)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

na_siteid_slope <- unique(slope[is.na(slope$linear_slope), ]$siteid)
sd_num <- analysis_dataset %>%
group_by(siteid) %>%
summarise(across(where(is.numeric), ~(sd(.x))))
sd_num %>%
filter(siteid %in% na_siteid_slope) %>%
summary(.)
#> siteid species_nb total_abundance log_total_abundance
#> Length:0 Min. : NA Min. : NA Min. : NA
#> Class :character 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Mode :character Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA
#> log_species_nb year latitude longitude chao_richness
#> Min. : NA Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA Max. : NA
#> chao_shannon chao_simpson chao_evenness hillebrand total
#> Min. : NA Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA Max. : NA
#> appearance disappearance jaccard horn chao
#> Min. : NA Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA Max. : NA
#> total_abundance_int shannon simpson inv_simpson evenness
#> Min. : NA Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA Max. : NA
#> jaccard_dis nestedness turnover dist_up_km ord_stra
#> Min. : NA Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA Max. : NA
#> dis_m3_pyr ele_mt_cav slp_dg_cav tmp_dc_cyr pac_pc_cse
#> Min. : NA Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA Max. : NA
#> urb_pc_cse hft_ix_c93 hft_ix_c09 riv_str_rc1 riv_str_rc2
#> Min. : NA Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA Max. : NA
#> first_year last_year span chao_richness_tps species_nb_tps
#> Min. : NA Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA Max. : NA
#> total_abundance_tps jaccard_scaled jaccard_dis_scaled nestedness_scaled
#> Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA
#> turnover_scaled hillebrand_scaled hillebrand_dis hillebrand_dis_scaled
#> Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA
#> appearance_scaled disappearance_scaled evenness_scaled
#> Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA
#> chao_richness_tps_scaled species_nb_tps_scaled scaled_dist_up_km
#> Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA
#> log_dist_up_km scaled_tmp_dc_cyr tmp_c_cyr scaled_tmp_c_cyr
#> Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA
#> hft_ix_c9309_percent log_hft_ix_c9309_percent hft_ix_c9309_ratio
#> Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA
#> hft_ix_c9309_log_ratio hft_ix_c9309_diff hft_ix_c9309_diff_scaled
#> Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA
#> one year_nb log1_year_nb total_abundance_scaled
#> Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA
- Tremendous trends in abundance:
high_abundance_slope <- slope_df %>%
arrange(desc(abs(total_abundance)))
high_abundance_slope
#> # A tibble: 5,516 x 22
#> siteid total_abundance log_total_abundance species_nb log_species_nb
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 S1661 -1726. -66.9 -0.130 -1.49
#> 2 S524 1027. 10.6 0.174 0.850
#> 3 S2779 941. 16.8 0.480 2.01
#> 4 S11735 -756. -0.971 0.245 1.31
#> 5 S2371 627. 11.0 0.408 1.48
#> 6 S6466 -472. -7.44 0.309 1.38
#> 7 S3655 436. 11.6 0.632 2.86
#> 8 S435 428. 31.4 0.655 4.56
#> 9 S1199 405. -16.3 -0.0588 -0.718
#> 10 S5180 398. 3.91 0.207 0.822
#> # ... with 5,506 more rows, and 17 more variables: chao_richness <dbl>,
#> # chao_shannon <dbl>, chao_simpson <dbl>, chao_evenness <dbl>, jaccard <dbl>,
#> # horn <dbl>, chao <dbl>, hillebrand <dbl>, total <dbl>, appearance <dbl>,
#> # disappearance <dbl>, evenness <dbl>, shannon <dbl>, simpson <dbl>,
#> # jaccard_dis <dbl>, nestedness <dbl>, turnover <dbl>
site_desc_slope_total_abundance <-
high_abundance_slope %>%
.[["siteid"]]
- The most tremendous temporal trends in abundance are located in
Illinois and France
ti <- filtered_dataset$location %>%
filter(siteid %in% site_desc_slope_total_abundance[1:20]) %>%
st_as_sf(coords = c("longitude", "latitude"), crs = 2154)
sp::plot(ne_countries())
#> Warning in wkt(obj): CRS object has no comment
sp::plot(st_geometry(ti), add = TRUE, col = "red", pch = 19)

filtered_dataset$abun_rich_op %>%
filter(siteid %in% site_desc_slope_total_abundance[1:20]) %>%
group_by(siteid) %>%
summarise(across(c("total_abundance", "species_nb"), list(median = median, mean = mean)))
#> # A tibble: 20 x 5
#> siteid total_abundance_med~ total_abundance~ species_nb_medi~ species_nb_mean
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 S1025 2534. 2960. 21 19.1
#> 2 S11735 5948 17416. 18 18.7
#> 3 S1199 0.443 5760. 8 8.67
#> 4 S1661 10700 13043. 9 9
#> 5 S1698 314. 866 14 12
#> 6 S1787 2790. 2498. 16.5 15.7
#> 7 S2371 6251 7982. 28 27.8
#> 8 S2779 2510 7132. 23 24.2
#> 9 S3655 2298 3689. 23 23.3
#> 10 S3863 1827 3353. 13 13.6
#> 11 S430 1987 3201. 12.5 12.5
#> 12 S435 2847 3087. 16 15.4
#> 13 S436 1936 2506. 19 20
#> 14 S4501 2111 1924. 16 14.9
#> 15 S5180 4335 8431. 23 23.3
#> 16 S524 27971 22155. 26 25.3
#> 17 S5812 4968 12216. 24.5 25.2
#> 18 S6466 3166 6465. 27 25.7
#> 19 S8158 3610 5140. 24.5 24.7
#> 20 S8577 4474 7435. 28 26.4
p_high_slope_abun <- map(site_desc_slope_total_abundance[1:20],
~plot_community_data(
analysis_dataset %>%
filter(siteid == .x) %>%
select(siteid, year, total_abundance),
y = "total_abundance", x = "year",
smoothing_method = "lm"))
plot_grid(plotlist = p_high_slope_abun, ncol = 3)

s11735_com <- analysis_dataset %>%
filter(siteid == site_desc_slope_total_abundance[1]) %>%
arrange(year) %>%
select(siteid, year, total_abundance)
s11735_com
#> # A tibble: 7 x 3
#> siteid year total_abundance
#> <chr> <dbl> <dbl>
#> 1 S1661 2003 18900
#> 2 S1661 2005 14500
#> 3 S1661 2007 10700
#> 4 S1661 2008 47200
#> 5 S1661 2012 0.129
#> 6 S1661 2013 0.154
#> 7 S1661 2018 0.066
plot_community_data(s11735_com, y = "total_abundance", x = "year",
smoothing_method = "lm")

mean(s11735_com$total_abundance) - sd(s11735_com$total_abundance)
#> [1] -3883.247
mean(s11735_com$total_abundance) + 5 * sd(s11735_com$total_abundance)
#> [1] 97673.67
s11735_pop <- filtered_dataset$measurement %>%
filter(siteid == site_desc_slope_total_abundance[1]) %>%
arrange(year) %>%
select(siteid, year, species, abundance)
s11735_pop %>%
group_by(species) %>%
summarise(abun = sum(abundance)) %>%
arrange(desc(abun))
#> # A tibble: 13 x 2
#> species abun
#> <chr> <dbl>
#> 1 Leuciscus leuciscus 26600.
#> 2 Rutilus rutilus 25700.
#> 3 Gobio gobio 13600.
#> 4 Perca fluviatilis 11200.
#> 5 Alburnus alburnus 4800.
#> 6 Anguilla anguilla 3700.
#> 7 Squalius cephalus 2500.
#> 8 Esox lucius 2000.
#> 9 Tinca tinca 400
#> 10 Gymnocephalus cernua 300.
#> 11 Barbus barbus 200
#> 12 Cyprinus carpio 200
#> 13 Abramis brama 100.
plot_temporal_population(s11735_pop) +
theme(legend.position = "bottom")
#> Note: Using an external vector in selections is ambiguous.
#> i Use `all_of(species_var)` instead of `species_var` to silence this message.
#> i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
#> This message is displayed once per session.
#> Note: Using an external vector in selections is ambiguous.
#> i Use `all_of(y_var)` instead of `y_var` to silence this message.
#> i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
#> This message is displayed once per session.
#> Note: Using an external vector in selections is ambiguous.
#> i Use `all_of(species)` instead of `species` to silence this message.
#> i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
#> This message is displayed once per session.
#> Warning: Removed 117 rows containing missing values (position_stack).

- Second highest temporal trends in abundance:
s6466_com <- analysis_dataset %>%
filter(siteid == site_desc_slope_total_abundance[2]) %>%
arrange(year) %>%
select(siteid, year, total_abundance)
s6466_com
#> # A tibble: 9 x 3
#> siteid year total_abundance
#> <chr> <dbl> <dbl>
#> 1 S524 1978 1861
#> 2 S524 1983 2739
#> 3 S524 1987 12032
#> 4 S524 1989 34800
#> 5 S524 1991 44608
#> 6 S524 1995 17982
#> 7 S524 1997 27971
#> 8 S524 2001 28139
#> 9 S524 2005 29259
plot_community_data(s6466_com, y = "total_abundance", x = "year",
smoothing_method = "lm")

s6466_pop <- filtered_dataset$measurement %>%
filter(siteid == site_desc_slope_total_abundance[2]) %>%
arrange(year) %>%
select(siteid, year, species, abundance)
s6466_pop %>%
group_by(species) %>%
summarise(abun = sum(abundance)) %>%
arrange(desc(abun))
#> # A tibble: 55 x 2
#> species abun
#> <chr> <dbl>
#> 1 Dorosoma cepedianum 93893
#> 2 Dorosoma petenense 68072
#> 3 Aplodinotus grunniens 14574
#> 4 Ictalurus furcatus 9962
#> 5 Ictalurus punctatus 5411
#> 6 Alosa chrysochloris 3645
#> 7 Notropis wickliffi 545
#> 8 Pylodictis olivaris 465
#> 9 Morone saxatilis 444
#> 10 Pomoxis annularis 429
#> # ... with 45 more rows
plot_temporal_population(s6466_pop) +
theme(legend.position = "bottom")
#> Warning: Removed 1045 rows containing missing values (position_stack).

Relationships among temporal trends
Globally
cor_slope <- slope_df %>%
select(-siteid) %>%
cor(., method = "spearman")
cor_slope
#> total_abundance log_total_abundance species_nb
#> total_abundance 1.0000000000 0.8246307994 0.31807533
#> log_total_abundance 0.8246307994 1.0000000000 0.40301195
#> species_nb 0.3180753260 0.4030119484 1.00000000
#> log_species_nb 0.2892004606 0.4052594645 0.93758602
#> chao_richness 0.0508961756 0.1394440304 0.80433662
#> chao_shannon -0.0279137551 0.0477362851 0.59176245
#> chao_simpson -0.0338305788 0.0191291401 0.48346752
#> chao_evenness -0.0444857958 0.0009310158 0.17645002
#> jaccard -0.0491817382 -0.0517691046 -0.08847793
#> horn -0.0456593373 -0.0471057146 -0.07711548
#> chao -0.0027001396 -0.0262558468 -0.01843873
#> hillebrand -0.0203384951 -0.0071825278 -0.06012651
#> total 0.0719557567 0.0788475302 0.16716036
#> appearance 0.2268847971 0.2929139747 0.69631041
#> disappearance -0.1710111291 -0.2509080868 -0.58273933
#> evenness -0.1304340666 -0.0971770501 0.19029318
#> shannon -0.0004716725 0.0813750746 0.60038484
#> simpson -0.0141998364 0.0516382845 0.48153130
#> jaccard_dis 0.0491429875 0.0517029275 0.08864262
#> nestedness 0.0345864709 0.0430393920 0.14269762
#> turnover 0.0132785000 -0.0069467758 -0.04408210
#> log_species_nb chao_richness chao_shannon chao_simpson
#> total_abundance 0.28920046 0.05089618 -0.02791376 -0.03383058
#> log_total_abundance 0.40525946 0.13944403 0.04773629 0.01912914
#> species_nb 0.93758602 0.80433662 0.59176245 0.48346752
#> log_species_nb 1.00000000 0.79445957 0.64663632 0.55175635
#> chao_richness 0.79445957 1.00000000 0.69826636 0.56653555
#> chao_shannon 0.64663632 0.69826636 1.00000000 0.94157972
#> chao_simpson 0.55175635 0.56653555 0.94157972 1.00000000
#> chao_evenness 0.24923947 0.19437940 0.65705554 0.69265391
#> jaccard -0.09704400 -0.09635087 -0.07884636 -0.07917356
#> horn -0.08126924 -0.08544403 -0.06734601 -0.07002592
#> chao -0.03482762 -0.02808783 -0.02116922 -0.02728646
#> hillebrand -0.05626267 -0.06552551 -0.08700625 -0.09673906
#> total 0.17660072 0.16318394 0.12898031 0.11747598
#> appearance 0.73091884 0.60465541 0.48732557 0.42063453
#> disappearance -0.62969158 -0.48910181 -0.40547117 -0.34222041
#> evenness 0.28946517 0.30909021 0.71480914 0.76973418
#> shannon 0.66069142 0.66738929 0.94205155 0.90503802
#> simpson 0.56369541 0.54659410 0.86919337 0.88186259
#> jaccard_dis 0.09705166 0.09639941 0.07884984 0.07917159
#> nestedness 0.14811936 0.14094149 0.11311663 0.09161961
#> turnover -0.06162618 -0.04540189 -0.04807683 -0.03371042
#> chao_evenness jaccard horn chao
#> total_abundance -0.0444857958 -0.049181738 -0.045659337 -0.002700140
#> log_total_abundance 0.0009310158 -0.051769105 -0.047105715 -0.026255847
#> species_nb 0.1764500235 -0.088477927 -0.077115480 -0.018438735
#> log_species_nb 0.2492394721 -0.097043997 -0.081269239 -0.034827625
#> chao_richness 0.1943794001 -0.096350866 -0.085444026 -0.028087826
#> chao_shannon 0.6570555375 -0.078846363 -0.067346005 -0.021169218
#> chao_simpson 0.6926539065 -0.079173560 -0.070025916 -0.027286459
#> chao_evenness 1.0000000000 -0.003071098 -0.001254299 -0.006306812
#> jaccard -0.0030710979 1.000000000 0.987437694 0.551615788
#> horn -0.0012542986 0.987437694 1.000000000 0.610948564
#> chao -0.0063068119 0.551615788 0.610948564 1.000000000
#> hillebrand -0.0356622276 0.524436363 0.518092095 0.198387953
#> total 0.0133645166 -0.959438300 -0.941968794 -0.488360340
#> appearance 0.1560546085 -0.645639868 -0.628584932 -0.295500150
#> disappearance -0.1607903320 -0.578334396 -0.586403546 -0.279086762
#> evenness 0.7290746559 -0.025293259 -0.020764263 -0.029613862
#> shannon 0.5893911297 -0.075280020 -0.063088213 -0.028773640
#> simpson 0.6081360040 -0.060030935 -0.049908681 -0.034585324
#> jaccard_dis 0.0031064970 -0.999918518 -0.987354545 -0.550443420
#> nestedness 0.0032391978 -0.464244804 -0.418463495 -0.350854086
#> turnover -0.0125261217 -0.509148654 -0.535620536 -0.137217564
#> hillebrand total appearance disappearance
#> total_abundance -0.020338495 0.07195576 0.22688480 -0.17101113
#> log_total_abundance -0.007182528 0.07884753 0.29291397 -0.25090809
#> species_nb -0.060126512 0.16716036 0.69631041 -0.58273933
#> log_species_nb -0.056262666 0.17660072 0.73091884 -0.62969158
#> chao_richness -0.065525508 0.16318394 0.60465541 -0.48910181
#> chao_shannon -0.087006248 0.12898031 0.48732557 -0.40547117
#> chao_simpson -0.096739058 0.11747598 0.42063453 -0.34222041
#> chao_evenness -0.035662228 0.01336452 0.15605461 -0.16079033
#> jaccard 0.524436363 -0.95943830 -0.64563987 -0.57833440
#> horn 0.518092095 -0.94196879 -0.62858493 -0.58640355
#> chao 0.198387953 -0.48836034 -0.29550015 -0.27908676
#> hillebrand 1.000000000 -0.50582982 -0.35709332 -0.32864160
#> total -0.505829820 1.00000000 0.72545884 0.52692942
#> appearance -0.357093319 0.72545884 1.00000000 -0.06905901
#> disappearance -0.328641599 0.52692942 -0.06905901 1.00000000
#> evenness -0.063799955 0.03861596 0.18872876 -0.18425549
#> shannon -0.085863584 0.12539124 0.48782752 -0.41692496
#> simpson -0.075107040 0.10005122 0.40842726 -0.36037753
#> jaccard_dis -0.524527998 0.95951672 0.64569401 0.57838092
#> nestedness -0.164497385 0.45337385 0.29607741 0.05875095
#> turnover -0.353682040 0.49179103 0.37125410 0.57019055
#> evenness shannon simpson jaccard_dis
#> total_abundance -0.13043407 -0.0004716725 -0.01419984 0.049142987
#> log_total_abundance -0.09717705 0.0813750746 0.05163828 0.051702928
#> species_nb 0.19029318 0.6003848409 0.48153130 0.088642617
#> log_species_nb 0.28946517 0.6606914197 0.56369541 0.097051655
#> chao_richness 0.30909021 0.6673892861 0.54659410 0.096399411
#> chao_shannon 0.71480914 0.9420515514 0.86919337 0.078849845
#> chao_simpson 0.76973418 0.9050380214 0.88186259 0.079171595
#> chao_evenness 0.72907466 0.5893911297 0.60813600 0.003106497
#> jaccard -0.02529326 -0.0752800203 -0.06003094 -0.999918518
#> horn -0.02076426 -0.0630882128 -0.04990868 -0.987354545
#> chao -0.02961386 -0.0287736400 -0.03458532 -0.550443420
#> hillebrand -0.06379995 -0.0858635844 -0.07510704 -0.524527998
#> total 0.03861596 0.1253912449 0.10005122 0.959516724
#> appearance 0.18872876 0.4878275209 0.40842726 0.645694012
#> disappearance -0.18425549 -0.4169249616 -0.36037753 0.578380917
#> evenness 1.00000000 0.7759287696 0.85878029 0.025227287
#> shannon 0.77592877 1.0000000000 0.96084041 0.075250808
#> simpson 0.85878029 0.9608404079 1.00000000 0.059989068
#> jaccard_dis 0.02522729 0.0752508078 0.05998907 1.000000000
#> nestedness 0.04733133 0.1203281052 0.10014131 0.464283551
#> turnover -0.04622306 -0.0635190179 -0.06433463 0.509189404
#> nestedness turnover
#> total_abundance 0.034586471 0.013278500
#> log_total_abundance 0.043039392 -0.006946776
#> species_nb 0.142697621 -0.044082096
#> log_species_nb 0.148119356 -0.061626179
#> chao_richness 0.140941491 -0.045401891
#> chao_shannon 0.113116630 -0.048076828
#> chao_simpson 0.091619608 -0.033710423
#> chao_evenness 0.003239198 -0.012526122
#> jaccard -0.464244804 -0.509148654
#> horn -0.418463495 -0.535620536
#> chao -0.350854086 -0.137217564
#> hillebrand -0.164497385 -0.353682040
#> total 0.453373854 0.491791029
#> appearance 0.296077412 0.371254096
#> disappearance 0.058750950 0.570190549
#> evenness 0.047331331 -0.046223059
#> shannon 0.120328105 -0.063519018
#> simpson 0.100141306 -0.064334629
#> jaccard_dis 0.464283551 0.509189404
#> nestedness 1.000000000 -0.389647237
#> turnover -0.389647237 1.000000000
png(here("doc", "fig", "p_cor_slope.png"))
corrplot::corrplot(cor_slope[c("chao_richness", "species_nb", "log_species_nb",
"total_abundance", "log_total_abundance"), ], type = "upper", tl.srt = 45, diag = FALSE)
dev.off()
#> pdf
#> 2
png(here("doc", "fig", "p_cor_slope_tot.png"))
corrplot::corrplot(cor_slope, type = "upper", tl.srt = 35, diag = FALSE)
dev.off()
#> pdf
#> 2
corrplot::corrplot(cor_slope[c("chao_richness", "species_nb", "log_species_nb",
"total_abundance", "log_total_abundance"), ], type = "upper", tl.srt = 45, diag = FALSE)

corrplot::corrplot(cor_slope, type = "upper", tl.srt = 35, diag = FALSE)

ti <- expand.grid(
resp1 = unique(slope$response),
resp2 = unique(slope$response)
) %>%
filter(resp2 != resp1) %>%
filter(
resp1 %in% c("chao_richness", "species_nb", "log_species_nb",
"total_abundance", "log_total_abundance")) %>%
mutate_all(as.character) %>%
arrange(resp1)
test <- map2(ti$resp1, ti$resp2,
function(x, y) {
bi <- slope %>%
filter(response %in% c(x, y)) %>%
select(siteid, response, linear_slope) %>%
pivot_wider(names_from = "response", values_from = "linear_slope")
return(bi)
}
)
p_trends_trends <- map(test, function(x) {
l <- colnames(x)
x %>%
ggplot(aes(x = !!sym(l[2]), y = !!sym(l[3]))) +
geom_point() +
geom_smooth(method = "loess") +
labs(
x = get_var_replacement()[l[2]],
y = get_var_replacement()[l[3]]
)
})
names(p_trends_trends) <- map_chr(test, ~colnames(.x)[2])
Total abundance
plot_grid(
plotlist = p_trends_trends[names(p_trends_trends) %in% "total_abundance"],
ncol = 3)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

Species richness
plot_grid(
plotlist = p_trends_trends[names(p_trends_trends) %in% "species_nb"],
ncol = 3)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

Log species richness
plot_grid(
plotlist = p_trends_trends[names(p_trends_trends) %in% "log_species_nb"],
ncol = 3
)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

- Trends of total abundance are linked to the one of species richness but not to
shannon and simpson, suggesting that the trends of abundance are more link to
trends of rare species. Logic: gain and loss of species should concern rare
species.
- Species richness decrease linked to disappearance, but not increase
- Species richness increase linked to appearance, but not decrease
Turnover vs nestedness:
test <- slope_df %>%
mutate(
nes_sup_tu = nestedness > turnover,
nes_inf_tu = nestedness < turnover,
rich_inc = species_nb >= 0
)
sum(test$nes_sup_tu)
#> [1] 2608
sum(test$nes_inf_tu)
#> [1] 2603
test %>%
tabyl(rich_inc, nes_sup_tu) %>%
adorn_totals("row") %>%
adorn_title(placement = "top")
#> nes_sup_tu
#> rich_inc FALSE TRUE
#> FALSE 1408 1037
#> TRUE 1500 1571
#> Total 2908 2608
median_site_rich <- analysis_dataset %>%
select(c("siteid", starts_with("chao_"), "total_abundance")) %>%
group_by(siteid) %>%
summarise(across(where(is.double), median)) %>%
rename_with(~paste0("med_", .x), where(is.double))
slope_df_med <- slope_df %>%
left_join(median_site_rich, by = "siteid") %>%
mutate(
nes_sup_tu = nestedness > turnover,
nes_inf_tu = nestedness < turnover,
rich_inc = species_nb >= 0
)
slope_df_med %>%
ggplot(aes(
x = log_species_nb,
y = jaccard_dis,
shape = nes_sup_tu,
color = nes_sup_tu
)) +
geom_point() +
geom_smooth(method = "gam")
#> `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'

p_baselga_jaccard <- map(c("nestedness", "turnover"),
~slope_df_med %>%
ggplot(aes(
x = jaccard_dis,
y = !!sym(.x),
color = log(med_chao_richness))) +
geom_point() +
scale_color_viridis() +
geom_abline(intercept = 0, slope = 1) +
geom_smooth(method = "loess")
)
p_baselga_jaccard[[3]] <-
slope_df_med %>%
ggplot(aes(
x = jaccard_dis,
y = nestedness - turnover,
color = log(med_chao_richness))) +
geom_point() +
scale_color_viridis() +
geom_abline(intercept = 0, slope = 1) +
geom_smooth(method = "loess")
plot_grid(plotlist = p_baselga_jaccard, ncol = 3)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

slope_df_med %>%
ggplot(aes(x = med_chao_richness, y = jaccard_dis)) +
geom_point()

Look at strong turnover trends
high_jaccard_slope <-
slope_df %>%
select(siteid, total_abundance, species_nb, jaccard) %>%
arrange(desc(jaccard))
high_jaccard_slope
#> # A tibble: 5,516 x 4
#> siteid total_abundance species_nb jaccard
#> <chr> <dbl> <dbl> <dbl>
#> 1 S6544 0.257 -0.0769 0.0385
#> 2 S8076 -1.99 -0.0383 0.0307
#> 3 S7264 1.16 0.0600 0.0300
#> 4 S6550 1.92 0.116 0.0286
#> 5 S5095 -3.37 -0.0857 0.0286
#> 6 S7806 0.760 -0.0394 0.0285
#> 7 S2742 1.08 -0.00909 0.0258
#> 8 S6783 0.508 -0.0500 0.0250
#> 9 S450 120. 0.870 0.0247
#> 10 S12595 13.6 -0.0421 0.0228
#> # ... with 5,506 more rows
site_desc_jaccard_slope <-
high_jaccard_slope %>%
.[["siteid"]]
ti <- filtered_dataset$location %>%
filter(siteid %in% site_desc_jaccard_slope[1:20]) %>%
st_as_sf(coords = c("longitude", "latitude"), crs = 2154)
sp::plot(ne_countries())
#> Warning in wkt(obj): CRS object has no comment
sp::plot(st_geometry(ti), add = TRUE, col = "red", pch = 19)

filtered_dataset$abun_rich_op %>%
filter(siteid %in% site_desc_jaccard_slope[1:20]) %>%
group_by(siteid) %>%
summarise(across(c("total_abundance", "species_nb"),
list(median = median, mean = mean))) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "bordered", "hover")) %>%
scroll_box(width = "100%", height = "400px")
|
siteid
|
total_abundance_median
|
total_abundance_mean
|
species_nb_median
|
species_nb_mean
|
|
S12595
|
347.000
|
341.62632
|
1.0
|
1.368421
|
|
S2742
|
41.900
|
57.09091
|
2.0
|
2.000000
|
|
S2831
|
170.500
|
180.15455
|
1.0
|
1.272727
|
|
S450
|
877.500
|
744.66667
|
12.5
|
9.166667
|
|
S5095
|
37.900
|
52.81667
|
1.0
|
1.333333
|
|
S5190
|
56.600
|
82.35714
|
2.0
|
1.857143
|
|
S5514
|
25.900
|
23.58000
|
1.5
|
1.700000
|
|
S5672
|
6.400
|
41.36000
|
2.0
|
1.800000
|
|
S6544
|
15.000
|
18.06154
|
1.0
|
1.307692
|
|
S6550
|
6.650
|
14.55714
|
2.0
|
2.214286
|
|
S6783
|
10.100
|
12.24667
|
1.0
|
1.200000
|
|
S6936
|
22.800
|
26.62500
|
1.0
|
1.333333
|
|
S7264
|
11.300
|
10.94000
|
1.0
|
1.400000
|
|
S7288
|
8.100
|
15.01667
|
3.0
|
3.333333
|
|
S7342
|
13.039
|
13.71067
|
4.0
|
3.833333
|
|
S7532
|
1.900
|
1.84000
|
1.0
|
1.400000
|
|
S7806
|
10.200
|
10.95000
|
1.0
|
1.222222
|
|
S8076
|
65.000
|
72.50000
|
2.0
|
2.250000
|
|
S8447
|
10.700
|
11.00000
|
4.0
|
3.750000
|
|
S998
|
34.100
|
36.39412
|
1.0
|
1.235294
|
p_high_jaccard_slope <- map(site_desc_jaccard_slope[1:20],
~plot_community_data(
analysis_dataset %>%
filter(siteid == .x) %>%
select(siteid, year, jaccard),
y = "jaccard", x = "year",
smoothing_method = "lm"))
plot_grid(plotlist = p_high_jaccard_slope, ncol = 3)

- Hum, do not understand why community are more similar at the end of the time-steps
p_pop_high_jaccard <-
map(site_desc_jaccard_slope[1:20],
~plot_temporal_population(
filtered_dataset$measurement %>%
filter(siteid == .x)
) +
theme(legend.position = "none")
)
plot_grid(plotlist = p_pop_high_jaccard, ncol = 3)
#> Warning: Removed 16 rows containing missing values (position_stack).
#> Warning: Removed 10 rows containing missing values (position_stack).
#> Warning: Removed 15 rows containing missing values (position_stack).
#> Warning: Removed 3 rows containing missing values (position_stack).
#> Warning: Removed 24 rows containing missing values (position_stack).
#> Warning: Removed 8 rows containing missing values (position_stack).
#> Warning: Removed 50 rows containing missing values (position_stack).
#> Warning: Removed 30 rows containing missing values (position_stack).
#> Warning: Removed 4 rows containing missing values (position_stack).
#> Warning: Removed 24 rows containing missing values (position_stack).
#> Warning: Removed 18 rows containing missing values (position_stack).
#> Warning: Removed 36 rows containing missing values (position_stack).
#> Warning: Removed 2 rows containing missing values (position_stack).
#> Warning: Removed 45 rows containing missing values (position_stack).

- Lets see communities that becomes more dissimilar:
p_jaccard_low_jaccard <-
map(site_desc_jaccard_slope[(length(site_desc_jaccard_slope)-20):length(site_desc_jaccard_slope)],
~plot_community_data(
analysis_dataset %>%
filter(siteid == .x) %>%
select(siteid, year, jaccard),
y = "jaccard", x = "year",
smoothing_method = "lm"))
plot_grid(plotlist = p_jaccard_low_jaccard, ncol = 3)

p_pop_low_jaccard <-
map(site_desc_jaccard_slope[
(length(site_desc_jaccard_slope)-20):length(site_desc_jaccard_slope)],
~plot_temporal_population(
filtered_dataset$measurement %>%
filter(siteid == .x)
) +
theme(legend.position = "none")
)
plot_grid(plotlist = p_pop_low_jaccard, ncol = 3)
#> Warning: Removed 132 rows containing missing values (position_stack).
#> Warning: Removed 15 rows containing missing values (position_stack).
#> Warning: Removed 20 rows containing missing values (position_stack).
#> Warning: Removed 25 rows containing missing values (position_stack).
#> Warning: Removed 13 rows containing missing values (position_stack).
#> Warning: Removed 18 rows containing missing values (position_stack).
#> Warning: Removed 9 rows containing missing values (position_stack).
#> Warning: Removed 12 rows containing missing values (position_stack).
#> Removed 12 rows containing missing values (position_stack).
#> Warning: Removed 78 rows containing missing values (position_stack).
#> Warning: Removed 18 rows containing missing values (position_stack).
#> Warning: Removed 10 rows containing missing values (position_stack).
#> Warning: Removed 102 rows containing missing values (position_stack).
#> Warning: Removed 12 rows containing missing values (position_stack).
#> Warning: Removed 45 rows containing missing values (position_stack).
#> Warning: Removed 100 rows containing missing values (position_stack).
#> Warning: Removed 25 rows containing missing values (position_stack).
#> Warning: Removed 24 rows containing missing values (position_stack).
#> Warning: Removed 30 rows containing missing values (position_stack).
#> Warning: Removed 105 rows containing missing values (position_stack).

Comparison rigal, lm
rigal_trends_df <- map2_dfr(
rigal_trends, names(rigal_trends),
~.x %>% mutate(variable = .y)
) %>%
select(variable, siteid, linear_slope) %>%
pivot_wider(names_from = "variable", values_from = "linear_slope") %>%
mutate(type = "rigal")
slope_comp <- list()
slope_comp$rigal <- rigal_trends_df[order(rigal_trends_df$siteid), ]
slope_comp$simple_lm <- slope_df[order(slope_df$siteid), ]
stopifnot(all(slope_comp$rigal$siteid == slope_comp$simple_lm$siteid))
slope_comp$diff <- tibble(
siteid = slope_comp$rigal$siteid)
for (i in var_temporal_trends) {
slope_comp$diff[[i]] <- slope_comp$rigal[[i]] - slope_comp$simple_lm[[i]]
}
old_par <- par()
par(mfrow = c(3, 2))
for (i in var_temporal_trends) {
plot(
slope_comp$rigal[[i]] ~
slope_comp$simple_lm[[i]],
ylab = "Rigal (polynomial degree 2)",
xlab = "Simple LM"
)
abline(0, 1)
title(i)
}



par(old_par)
#> Warning in par(old_par): graphical parameter "cin" cannot be set
#> Warning in par(old_par): graphical parameter "cra" cannot be set
#> Warning in par(old_par): graphical parameter "csi" cannot be set
#> Warning in par(old_par): graphical parameter "cxy" cannot be set
#> Warning in par(old_par): graphical parameter "din" cannot be set
#> Warning in par(old_par): graphical parameter "page" cannot be set

# Compare Rigal and simple lm for few sites
plot_rigal_raw_data <- function(
site = NULL,
response = NULL,
dataset = NULL,
rigal_trends = NULL
) {
raw_data <- dataset %>%
filter(siteid == site)
rigal_resp <- rigal_trends[[response]]
rigal_resp_site <- rigal_resp[rigal_resp$siteid == site, ]
plot(raw_data[[response]]~raw_data[["year"]], ylab = response, xlab = "Year")
abline(rigal_resp_site$intercept, rigal_resp_site$linear_slope, col
= "red")
abline(lm(raw_data[[response]]~raw_data[["year"]]), col = "green")
title(main = site)
}
plot_rigal_raw_data(site = "S2672", response = "log_total_abundance", dataset =
analysis_dataset, rigal_trends = rigal_trends)

slope_comp$diff <- slope_comp$diff %>%
arrange(desc(abs(log_total_abundance)))
old_par <- par()
par(mfrow = c(3, 2))
for (i in slope_comp$diff[1:6,]$siteid) {
plot_rigal_raw_data(
site = i,
response = "log_total_abundance",
dataset = analysis_dataset,
rigal_trends = rigal_trends
)
}

par(old_par)
#> Warning in par(old_par): graphical parameter "cin" cannot be set
#> Warning in par(old_par): graphical parameter "cra" cannot be set
#> Warning in par(old_par): graphical parameter "csi" cannot be set
#> Warning in par(old_par): graphical parameter "cxy" cannot be set
#> Warning in par(old_par): graphical parameter "din" cannot be set
#> Warning in par(old_par): graphical parameter "page" cannot be set
Classification of temporal trends
classification_station <- map2(rigal_trends, names(rigal_trends),
~.x %>%
ggplot(aes(x = shape_class)) +
geom_bar() +
labs(x = "Trajectory classification", y = "Nb station",
title = get_var_replacement()[.y]) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
)
plot_grid(plotlist = classification_station)

p_classif_abun_rich <- plot_grid(
plotlist = classification_station[
names(classification_station)[str_detect(names(classification_station),
"abundance|richness|species")]])
save_plot(
here("doc", "fig", "classif_trends_abun_rich.png"),
p_classif_abun_rich
)
p_classif_hill_evenness <- plot_grid(
plotlist = classification_station[
names(classification_station)[str_detect(names(classification_station),
"shannon|simpson|evennes")]])
save_plot(
here("doc", "fig", "classif_trends_hill_evenness.png"),
p_classif_hill_evenness
)
p_classif_turnover <- plot_grid(
plotlist = classification_station[
names(classification_station)[names(classification_station) %in%
c("jaccard", "horn", "chao", "hillebrand", "total", "appearance",
"disappearance")]])
save_plot(
here("doc", "fig", "classif_trends_classif_hill_evenness.png"),
p_classif_turnover
)
ti <- map_dfr(rigal_trends, ~tabyl_df(x = .x, group = "direction"),
.id = "response"
)
ti %>%
filter(direction != "Total") %>%
select(response, direction, percent) %>%
mutate(response = str_replace_all(response, get_var_replacement())) %>%
pivot_wider(names_from = "direction", values_from = "percent") %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "bordered", "hover"))
|
response
|
stable
|
increase
|
decrease
|
|
Total abundance
|
82.7%
|
10.0%
|
7.3%
|
|
Log Total Turnover (jaccard) (codyn) abundance
|
80.6%
|
11.8%
|
7.6%
|
|
Species richness
|
81.1%
|
11.5%
|
5.7%
|
|
Log species richness
|
79.5%
|
11.1%
|
5.2%
|
|
Chao species richness
|
83.2%
|
9.4%
|
6.1%
|
|
Chao Shannon
|
80.7%
|
9.8%
|
5.7%
|
|
Chao Simpson
|
81.7%
|
8.8%
|
5.8%
|
|
Chao Evenness
|
85.1%
|
6.3%
|
4.9%
|
|
Jaccard (binary, similarity)
|
70.4%
|
0.6%
|
27.2%
|
|
Horn (binary, similarity)
|
71.6%
|
0.6%
|
26.0%
|
|
Chao (binary, similarity)
|
72.4%
|
2.6%
|
12.6%
|
|
SER_a (rel abundance)
|
73.5%
|
0.7%
|
24.5%
|
|
Total Turnover (jaccard) (codyn)
|
68.2%
|
26.0%
|
0.3%
|
|
Appearance
|
68.3%
|
19.7%
|
0.4%
|
|
disAppearance
|
61.2%
|
14.8%
|
0.2%
|
|
Evenness
|
82.5%
|
7.6%
|
6.2%
|
|
Shannon
|
80.0%
|
10.2%
|
6.1%
|
|
Simpson
|
81.5%
|
9.0%
|
5.8%
|
|
Jaccard (binary, dissimilarity)
|
67.4%
|
26.8%
|
0.3%
|
|
Nestedness (jaccard)
|
80.3%
|
13.1%
|
0.5%
|
|
Turnover (jaccard)
|
53.9%
|
12.2%
|
0.2%
|
#3.3 Log-linear model:logYi=α+βXi+εiIn
#the log-linear model, the literal interpretation of the estimated
#coefficientˆβis that a one-unitincrease inXwill produce an expected increase in
#logYofˆβunits.
#In terms ofYitself, this
#meansthat the expected value ofYis multiplied byeˆβ. So in terms of
#effects of changes inXonY(unlogged):•Each 1-unit increase inXmultiplies the
#expected value ofYbyeˆβ.•To compute the effects onYof another change inXthan an
#increase of one unit, call thischangec, we need to includecin the exponent. The
#effect of ac-unit increase inXis tomultiply the expected value ofYbyecˆβ. So the
#effect for a 5-unit increase inXwould bee5ˆβ.•For small values ofˆβ,
#approximatelyeˆβ≈1+ˆβ. We can use this for the following approxima-tion for a
#quick interpretation of the coefficients: 100·ˆβis the expected percentage
#changeinYfor a unit increase inX. For instance forˆβ=.06,e.06≈1.06, so a 1-unit
#change inXcorresponds to (approximately) an expected increase inYof 6%
#https://data.library.virginia.edu/interpreting-log-transformations-in-a-linear-model/
hist_linear_plot <- map2(rigal_trends, names(rigal_trends),
~.x %>%
ggplot(aes(x = linear_slope)) +
geom_histogram() +
labs(x = "Linear slope", y = "Number of station",
title = get_var_replacement()[.y]) +
ylim(c(0, 1000))
)
plot_grid(plotlist = hist_linear_plot)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 3 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 2 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 2 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 2 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 2 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 2 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 2 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).

type_response <- list(
abun_rich = str_detect(names(classification_station),
"abundance|richness|species"),
hill_evenness = str_detect(names(classification_station),
"shannon|simpson|evennes"),
turnover = names(classification_station) %in%
c("jaccard", "horn", "chao", "hillebrand", "total", "appearance",
"disappearance")
)
map2(type_response, names(type_response),
function(x, y) {
p <- plot_grid(
plotlist = hist_linear_plot[
names(hist_linear_plot)[x]
]
)
save_plot(
here("doc", "fig", paste0("hist_linear_slope_", y, ".png")),
p
)
}
)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 3 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 2 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 2 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 2 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 2 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 2 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> $abun_rich
#> [1] "L:/alain/RivFishTimeBiodiversityFacets/doc/fig/hist_linear_slope_abun_rich.png"
#>
#> $hill_evenness
#> [1] "L:/alain/RivFishTimeBiodiversityFacets/doc/fig/hist_linear_slope_hill_evenness.png"
#>
#> $turnover
#> [1] "L:/alain/RivFishTimeBiodiversityFacets/doc/fig/hist_linear_slope_turnover.png"
Spatialization of temporal trends
rigal_trends_df_loc <-
rigal_trends_df %>%
left_join(
select(filtered_dataset$location, siteid, ecoregion, country),
by = c("siteid"))
- No obvious geographic patterns
rigal_trends_df_loc %>%
ggplot(aes(x = country, y = log_species_nb, fill = ecoregion)) +
geom_boxplot()

rigal_trends_df_loc %>%
ggplot(aes(x = country, y = log_total_abundance, fill = ecoregion)) +
geom_boxplot()

rigal_trends_df_loc %>%
ggplot(aes(x = country, y = total, fill = ecoregion)) +
geom_boxplot()

rigal_trends_df_loc %>%
ggplot(aes(x = country, y = jaccard, fill = ecoregion)) +
geom_boxplot()

rigal_trends_df_loc %>%
ggplot(aes(x = country, y = jaccard_dis, fill = ecoregion)) +
geom_boxplot()

rigal_trends_df_loc %>%
ggplot(aes(x = country, y = hillebrand, fill = ecoregion)) +
geom_boxplot()

rigal_trends_df_loc <- filtered_dataset$location %>%
left_join(rigal_trends_df, by = "siteid") %>%
st_as_sf(coords = c("longitude", "latitude"),
crs = 4326)
library(colorspace)
div_pal <- colorspace::sequential_hcl(11, "plasma")
world <- ne_countries(returnclass = "sf")
bb <- st_bbox(rigal_trends_df_loc)
map_world_trends(df = rigal_trends_df_loc, y = "log_species_nb")

map_world_trends(df = rigal_trends_df_loc, y = "jaccard")

map_world_trends(df = rigal_trends_df_loc, y = "log_total_abundance")

rigal_trends_df %>%
summarise(across(where(is.numeric), median, na.rm = T)) %>%
unlist()
#> total_abundance log_total_abundance species_nb log_species_nb
#> 9.553162e-02 5.582744e-01 3.893853e-03 1.305570e-01
#> chao_richness chao_shannon chao_simpson chao_evenness
#> 2.787916e-03 3.308454e-03 2.601314e-03 0.000000e+00
#> jaccard horn chao hillebrand
#> -1.527040e-02 -1.041153e-02 -3.884971e-17 -1.343464e-02
#> total appearance disappearance evenness
#> 1.195970e-02 5.364436e-03 3.065584e-03 2.796551e-05
#> shannon simpson jaccard_dis nestedness
#> 8.848320e-04 2.939331e-04 1.527040e-02 5.342415e-03
#> turnover
#> 2.203876e-03
median(exp(rigal_trends_df$log_total_abundance) - 1) * 100
#> [1] 74.76598
sd(exp(rigal_trends_df$log_total_abundance) - 1) * 100
#> [1] 1.494478e+35
median(exp(na.omit(rigal_trends_df$log_species_nb)) - 1) * 100
#> [1] 13.94629
sd(exp(na.omit(rigal_trends_df$log_species_nb)) - 1) * 100
#> [1] 1.551331e+19
Sensibility analysis of temporal trends
slope_avg3y <- map_dfr(rigal_trends_avg3y[
! names(rigal_trends_avg3y) %in% c("jaccard_dis", "nestedness", "turnover")], ~.x %>%
select(siteid, linear_slope),
.id = "response"
)
- Number of sites that are in the different datasets:
u_site <- map(list(classic = slope, avg_3y = slope_avg3y), ~unique(.x$siteid))
sum(u_site[[1]] %in% u_site[[2]])
#> [1] 5516
sum(u_site[[2]] %in% u_site[[1]])
#> [1] 5516
comparison_baseline <- rbind(
slope %>%
mutate(type = "classic") %>%
filter(siteid %in% u_site[["avg_3y"]]),
slope_avg3y %>%
mutate(type = "avg_3y")
) %>%
pivot_wider(names_from = type, values_from = linear_slope)
- It is not too bad, but we can see some differences:
- slope negatively biaised for
log_species_nb in
3 years baseline
- Bunch of trends in appearance/disappearance trends disappears with
3 years baseline
- But note that in three years baseline, there are three points less data
comparison_baseline %>%
ggplot(aes(x = classic, y = avg_3y)) +
geom_point() +
geom_abline(intercept = 0, slope = 1) +
facet_wrap(~response, scale = "free")
#> Warning: Removed 35396 rows containing missing values (geom_point).

Environment
col_to_keep <- get_river_atlas_significant_var()
names(col_to_keep) <- NULL
riv <- riveratlas_site %>%
clean_names() %>%
select(all_of(c("siteid", col_to_keep))) %>%
st_drop_geometry()
slp_riv <- slope_df %>%
select(siteid, log_species_nb, evenness, turnover, nestedness, hillebrand,
appearance, disappearance, jaccard) %>%
left_join(riv, by = "siteid") %>%
left_join(select(filtered_dataset$location, siteid, ecoregion, main_bas, latitude, longitude), by = "siteid")
length(unique(slp_riv$main_bas))
#> [1] 309
slp_riv %>%
ggplot(aes(y = jaccard, x = as.factor(ord_stra), color = ecoregion)) +
geom_boxplot()

slp_riv %>%
ggplot(aes(y = jaccard, x = log(dist_up_km), color = ecoregion)) +
geom_point() +
geom_smooth(method = "lm")
#> `geom_smooth()` using formula 'y ~ x'
#> Warning: Removed 600 rows containing non-finite values (stat_smooth).
#> Warning: Removed 600 rows containing missing values (geom_point).

# slp_riv %>%
# ggplot(aes(y = jaccard, x = as.factor(ord_flow), color = ecoregion)) +
# geom_boxplot()
slp_riv %>%
ggplot(aes(y = jaccard, x = tmp_dc_cyr /10, color = ecoregion)) +
geom_point() +
geom_smooth(method = "lm")
#> `geom_smooth()` using formula 'y ~ x'
#> Warning: Removed 600 rows containing non-finite values (stat_smooth).
#> Removed 600 rows containing missing values (geom_point).

slp_riv %>%
ggplot(aes(y = hillebrand, x = as.factor(ord_stra), color = ecoregion)) +
geom_boxplot()

slp_riv %>%
ggplot(aes(y = hillebrand, x = log(dist_up_km), color = ecoregion)) +
geom_point() +
geom_smooth(method = "lm")
#> `geom_smooth()` using formula 'y ~ x'
#> Warning: Removed 600 rows containing non-finite values (stat_smooth).
#> Warning: Removed 600 rows containing missing values (geom_point).

slp_riv %>%
ggplot(aes(y = hillebrand, x = tmp_dc_cyr / 10, color = ecoregion)) +
geom_point() +
geom_smooth(method = "lm")
#> `geom_smooth()` using formula 'y ~ x'
#> Warning: Removed 600 rows containing non-finite values (stat_smooth).
#> Removed 600 rows containing missing values (geom_point).

slp_riv %>%
ggplot(aes(y = hillebrand, x = log(dist_up_km), color = as.factor(main_bas))) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
theme(legend.position = "none")
#> `geom_smooth()` using formula 'y ~ x'
#> Warning: Removed 600 rows containing non-finite values (stat_smooth).
#> Warning: Removed 600 rows containing missing values (geom_point).

slp_riv %>%
ggplot(aes(y = hillebrand, x = ord_stra, color = as.factor(main_bas))) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
theme(legend.position = "none")
#> `geom_smooth()` using formula 'y ~ x'
#> Warning: Removed 600 rows containing non-finite values (stat_smooth).
#> Removed 600 rows containing missing values (geom_point).

ti <- slp_riv %>%
nest_by(main_bas) %>%
filter(nrow(data) > 10 & sum(!is.na(data$dist_up_km)) > 10)
tu <- ti %>%
mutate(
m = list(lm(jaccard ~ log(dist_up_km), data = data)),
coef = list(broom::tidy(m))
) %>%
unnest(coef) %>%
filter(term == "log(dist_up_km)")
tu %>%
ggplot(aes(x = estimate)) +
geom_histogram()
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(lm(jaccard ~ log(dist_up_km), slp_riv))
#>
#> Call:
#> lm(formula = jaccard ~ log(dist_up_km), data = slp_riv)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -0.100741 -0.008875 0.002875 0.012125 0.057293
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -0.0139322 0.0007663 -18.181 < 2e-16 ***
#> log(dist_up_km) -0.0013307 0.0002148 -6.196 6.25e-10 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.0171 on 4914 degrees of freedom
#> (600 observations deleted due to missingness)
#> Multiple R-squared: 0.007753, Adjusted R-squared: 0.007551
#> F-statistic: 38.39 on 1 and 4914 DF, p-value: 6.25e-10
SpaMM model
knitr::opts_chunk$set(eval = FALSE)
tu <- ti %>%
mutate(
m = list(lm(appearance ~ ord_stra, data = data))
,
coef = list(broom.mixed::tidy(m))
) %>%
unnest(coef) %>%
filter(term == "ord_stra") %>%
ungroup()
stat_data <- slp_riv %>%
mutate(
appearance = scale(appearance),
ord_stra = as.factor(ord_stra),
main_bas = as.factor(main_bas),
log_dist_up_km = log(dist_up_km)
)
m0 <- spaMM::fitme(appearance ~ ord_stra, stat_data)
plot(m0)
m1 <- spaMM::fitme(appearance ~ ord_stra + (1 | main_bas),
stat_data)
m2 <- spaMM::fitme(
appearance ~ ord_stra + (1 + ord_stra | main_bas),
stat_data
)
m2d <- spaMM::fitme(
appearance ~ dist_up_km + (1 + dist_up_km | main_bas),
stat_data
)
summary(m2d)
m2d <- spaMM::fitme(
appearance ~ log_dist_up_km + (1 + log_dist_up_km | ecoregion/main_bas/siteid),
stat_data
)
summary(m2d)
plot(m2d)
tm <- stat_data %>%
na.omit() %>%
mutate(
pred_appearance = predict(m2d,
re.form = ~log_dist_up_km + (1 + log_dist_up_km | ecoregion))[,1]
)
tm %>%
ggplot(aes(x = log_dist_up_km, y = pred_appearance, color = ecoregion)) +
geom_line()
formula <- paste0(
"jaccard",
" ~ ",
"year", "+", "year:ord_stra", " + ", "(1 + year:ord_stra| ecoregion/main_bas/siteid)"
)
spaMM::fitme(formula = formula, data = analysis_dataset)
formula <- paste0(
"jaccard",
" ~ ",
"dis_up_km + tmp_dc_cyr + (1 + dist_up_km + tmp_dc_cyr | ecoregion/main_bas/siteid)"
)
m3 <- spaMM::fitme(scale(appearance) ~ ord_stra + (1 | main_bas) +
Matern(1 | longitude + latitude),
na.omit(slp_riv), control.HLfit = list(algebra = "decorr"))
tu <- ti %>%
mutate(
m = list(lm(disappearance ~ log(dist_up_km), data = data)),
coef = list(broom::tidy(m))
) %>%
unnest(coef) %>%
filter(term == "log(dist_up_km)") %>%
ungroup()
tu %>%
ggplot(aes(x = estimate)) +
geom_histogram()
summary(lm(disappearance ~ log(dist_up_km), slp_riv))
cor_slp_env <-
slp_riv[, !colnames(slp_riv) %in% c("siteid", "ecoregion", "main_bas")] %>%
na.omit() %>%
cor(., method = "spearman")
corrplot::corrplot(
cor_slp_env
)
rigal_trends_df %>%
ggplot(aes(x = jaccard , y = hillebrand)) +
geom_point() +
# geom_smooth(method = "gam") +
labs(x = "Jaccard trends (similarity, binary)", y = "Hillebrand trends
(similarity, relative abundance)")
library(INLA)
a3 <- inla(scale(appearance) ~ ord_stra + f(main_bas, model = "iid"),
family = "gaussian",
control.predictor = list(link = 1, compute = TRUE),
control.compute = list(cpo = TRUE, dic = TRUE, config = TRUE),
data = slp_riv)
PCA on temporal trends
slope_df <- rigal_slp_df
res.pca <- dudi.pca(select(slope_df, -siteid),
scannf = FALSE, # Hide scree plot
nf = 5 # Number of components kept in the results
)
fviz_eig(res.pca)
p_pca <- map(list(c(1,2), c(2, 3), c(1, 3)),
~fviz_pca_var(res.pca,
axes = .x,
col.var = "contrib", # Color by contributions to the PC
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Avoid text overlapping
) +
theme(legend.position = "none") + labs(title = "")
)
plot_grid(plotlist = p_pca, ncol = 1)
Reproducibility
Reproducibility receipt
## datetime
Sys.time()
## repository
if(requireNamespace('git2r', quietly = TRUE)) {
git2r::repository()
} else {
c(
system2("git", args = c("log", "--name-status", "-1"), stdout = TRUE),
system2("git", args = c("remote", "-v"), stdout = TRUE)
)
}
## session info
sessionInfo()